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1 Introduction 

The interest to electronic processes in disordered systems was greatly inspired by the fascinating 
disorder effects in semiconductors, including doped and amorphous ones [1], and in mesoscopic 
metallic systems [2]. Unlike the Bloch states in fully periodic systems, the electronic spectrum 
of disordered systems generally includes both extended and localized states, their coexistence 
being related to the competition between kinetic and potential energy of Fermi particles. The 
main consequence of this competition is the possibility for Anderson transition from metallic 
to insulating state at zero temperature and sufficiently strong disorder [3]. The best studied 
situation is that of non-interacting electrons (that is, the single-electron approximation) in a 
certain random field. Efficiency of single-electron theories for each type of electronic states 
in disordered systems (at the Fermi level ep, they are mostly extended in metals and mostly 
localized in semiconductors) is assured by the presence of a certain small parameter, such 
as Xp/£ (where Xf is the Fermi wavelength, £ the mean free path) in metals [2] and na 3 in 
semiconductors [4] (n is the concentration of charge carriers, a the lattice parameter). 

The following discussion is addressed to the doped semiconducting systems. In traditional 
semiconductors, typical values of na 3 do not exceed 10~ 4 -r- 1CT 6 . This is determined by a very 
great localization radius r ~ (20-h50)a of a single localized shallow state and the Mott criterion 
for metallization in the impurity band: nrjj ~ 0.02 [1]. At doping levels below this value, the 
single-electron approach yeilds in a finite Fermi density of (localized) states vp = z/(e^) and 
.£h ! m the Mott law for hopping conductivity vs temperature T: cr(T) oc {exp}(— BT~ 1 / 4 ) , B ps 
2.1(7'qZ/f)- 1 / 4 [5]. However, it was shown by Efros and Shklovskii [6] that account of Coulomb 
interactions between such shallow states leads to formation of a "soft gap" in the unperturbed 
density of states v F near the Fermi level so that: vie — ep) 2 /e 6 until \e — ep\ ~ A = e 3 Up 2 / k 3 ^ 2 
(where e is the electron charge and k, the static dielectric constant). Consequently, the Mott 
law is changed to: cr(T) oc exp(_B'T -1 / 2 ), B' « e/(/tr ) 1 / 2 at sufficiently low temperatures, 
T < T c w e^roVp/K 2 . The above theoretic dependencies are in a good agreement with the bulk 
of experimental data in traditional semiconductors. 
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A special class of doped materials, displaying semiconducting, metallic, superconducting 
and various magnetic phases, is comprised by the doped perovskite systems [7,8]. From the 
point of view of standard theory of semiconductors, these materials exhibit extremely high val- 
ues of na 3 ~ 0.1 4- 0.5, that is more than three orders of magnitude higher than those observed 
in traditional semiconductors. So the experimentally observed hopping type of conductance, 
at sufficiently high temperatures, in perovskite manganites Lni-^A^MnOs (where Ln is a lan- 
thanide, A an alkali-earth metal) with doping levels x ~ 0.3 [9] is quite a striking fact. It 
indicates that Fermi states in such materials remain localized even at so heavy doping and, 
from the before cited Mott criterion, the upper limit for localization radius should be estimated 
as r ~ 0.4a. This is an opposite limit to the traditional case of shallow dopants, hence a 
different evolution of the excitation spectrum can be expected. In particular, the effects of 
electron-electron interaction can be much more pronounced. 

Of course, real perovskite manganites possess many other peculiar properties, as spin- 
dependent kinetic energy (it is just this dependence that suppresses kinetic energy in the para- 
magnetic phase) and the related Zener mechanism of double exchange [10], strong coupling of 
charge carriers to Jahn-Teller deformations [11], possible formation of small spin polarons [12], 
charge localization [13] and charge ordering [14], etc. However all the above mechanisms are 
usually discussed within uniform (that is, completely ordered) models, while our main focus 
now is on the effects specific for extremely strong disorder. To this end, we propose a model 
approximation, starting from a set of strictly localized single-site electronic states, with ran- 
dom energy in each site formed by the superposition of classic Coulomb potentials from other 
(occupied) sites and from the fixed charged dopants (so that the overall electroneutrality is as- 
sured). At the next step, the kinetic energy is considered as a small perturbation, triggering the 
hops between the nearest neighbor occupied and empty sites (if the necessary energy difference 
is compensated by phonons, a.c. electric fields, etc.). This model is interesting firstly as an 
opposite limit to the usual situation when the disorder is treated as a small perturbation of an 
initially uniform system, and secondly as a new realization of strongly correlated many-body 
system. 

Below in Sec. 2 we define the model parameters and describe the numeric processes to 
seek the ground state and analyze various types of excitation spectra at zero temperature, as 
functions of system size L and doping concentration x. Then the size, shape and topology 
independent behavior of the spectra is rapidly attained with growing L (already at L = 10 -j- 12 
cell units), and the main findings are: 

(a) a non-ergodicity of the full phase space with a definite distribution of local energy minima 

and practically identical excitation spectra with respect to all typical local minima; 

(b) asymmetric deviations from the mean-field oc (e — e F ) 2 behavior for the density of single- 

particle excitations; 

(c) vanishing of the density of pair excitations in the limit of low excitation energy; 

(d) finiteness of the total density of many-body excitations in the same limit. 

The further development of the model in Sec. 3 involves finite temperatures and hopping dy- 
namics through the usual mechanism of deformation potential for electron-phonon coupling. 
Then the non-ergodic structure of the system phase space at zero temperature leads to its 
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Figure 1: a) Elementary cubic cell corresponding to the LaMnC>3 structure, b) A 2 x 2 x 3 
parallelepiped sample of cubic lattice with a random distribution of charged dopants and charge 
carriers at concentration x = 1/4. 

non-Markovian kinetics at finite temperatures. A special numeric algorithm is developed, 
simulating this statistical process and giving the temperature behavior of electronic specific 
heat, relaxation times and diffusion coefficient (the latter being related to the conductivity via 
the Einstein relation). In particular, for the typical concentration value x = 1/3 we found: 
<t(T) oc exp(— bT~ a ), with a ~ 0.87, different from the Efros and Shklovskii value 1/2. The 
conclusions and perspectives for the further studies of the present model are discussed in Sec. 
4. 

2 Description of model. Zero Temperature 

The crystalline structure of lanthanum manganite LaMnC>3 is close to ideal perovskite with 
the (paramagnetic) lattice parameter a ~ 3.9A(Fig. la). Substitution of trivalent La 3+ by 
divalent alcali-earth ion (say, Ca 2+ ) brings an extra hole to the system which resides on a 
nearby manganese site, changing its state from Mn 3+ to Mn 4+ . For a single dopant in the 
lattice, an arbitrarily weak tunneling will be sufficient to produce the hole state equally shared 
between eight Mn sites, nearest neighbors to the dopant Ca. However, for a finite doping, there 
appear random energy differences between these sites, and if these differences are greater than 
the tunneling amplitude (the kinetic energy), the hole will mainly occupy the lowest energy site. 
The localization is also favoured here by the suppression of kinetic energy: firstly due to the 
presence of intercalating oxygens between manganese sites and secondly due to the incoherence 
of manganese spins at higher temperatures. 

Referring to this situation we consider the model where the dopant ions occupy randomly 
the central sites in the simple cubic lattice with probability x < 1 (Fig. lb). Each dopant 
releases one charge carrier into the crystal and thus acquires a unit charge e of opposite sign. 
In neglectance of hops between lattice sites, each carrier occupy single lattice site and the total 
electronic energy includes only Coulomb contributions: 




(1) 



where U(n) = ^ n ,^ n c(n)u(|n - n'|), l^(n') = J2n'^n c(n')u(|n - n' - 5\), u{r) = e 2 /Kr, c(n) is 
an occupation number for charge carrier in the lattice site n = a(ni,n 2 ,n 3 ) (with integer rii), 
while occupation of the dopant site n + 5, with 5 = a(|, \, |), in the same cell is defined by 
d(n). In what follows we consider fixed random configuration d(n) of "frozen" dopants, then 
the system ground state corresponds to such adjustment of the configuration c(n) of carriers 
that E c is a minimum. 

For arbitrary configuration c(n), when a carrier is taken off from the site n (if occupied) or 
put into this site (if empty), the full energy respectively decreases or increases by e(n) = U(n) — 
V(n), which can be thus associated with single-particle excitations of the Fermi liquid theory. 
Here both U(n) and V(n) are random, but only V(n), determined by the fixed configuration of 
dopants, can be considered an usual local random field of Anderson's model [3] whereas U(n), 
determined by the variable configuration of carriers, is substantially non-local. Hence each e(n) 
essentially depends on the positions of all other carriers, and the total E c strongly differs from 
the sum of all e(n). In this situation, there is no evidence for unique energy minimum and the 
structure of phase space can be very complicate. However, the consideration is simplified if 
one takes in mind that any two configurations in this space, satisfying the same normalization 
condition, Eq. (2), can be connected by a sequence (a phase trajectory) of single-particle moves 
between nearest neighbor occupied and empty sites. In fact, the following treatment is limited 
just to this class of trajectories. 

Besides the above indicated single-particle energies e(n), the excitation spectrum includes 
also the so-called pair energies [6]: 

e(n,n / ) = e(n / )-e(n)- M '(|n-n'|), 

that is the energy change at moving a carrier from the occupied site n to empty site n'. 
The last, "excitonic" term in this expression, accounting for the correlation between different 
single-site energies, is just responsible for opening of the Coulomb gap in the single-particle 
spectrum z/( s _ p )(e(n)). If the correlations between different pair energies e(n, n') and e(n', n") 
are neglected, which corresponds to the mean-field approximation, a conclusion can be drawn 
that the pair spectrum v p (e(n, n')) is ungapped [6]. But, as will be seen from the exact numeric 
analysis below, in fact the density of such excitations also vanishes at e — > 0, and this situation 
may be supposed to exist for higher order excitations as well. 

In our numeric procedure we consider lattice samples in the form of finite parallelepipeds 
and apply the following algorithm. The initial configurations do({n}) and co({n}) are defined 
by assigning them independently random values or 1 for each site n, so that the normalization 
condition holds: 

Li L 2 L 3 L1+IL2+IL3 + I 

Y^J^J^d(n 1 + -,n 2 + -,n 3 + -)= ^ cfa^m) = N (2) 

rai=l 712=1 n3 = i ni=l n2=i 713=1 

related to the doping level x through N = [xLiL 2 L 3 ]. Then we choose, from all the nearest 
neighbor pairs of occupied sites n and empty sites n' , the pair n and n' corresponding to 
the minimum value of pair energy: e^ 9 ^ = e(n ,rio). If e^™ 9 ^ ^ s negative, we change 
from the configuration Co({n}) to a new configuration ci({n}), moving the carrier from n 
to n' Q . This process of single-particle moves is repeated m times, until we come to such a 
configuration c m ({n}) that ^in^ * s already positive. Then we search for the minimum 
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Figure 2: a) Single-particle density of states at T = and doping level x = 1/3. Insert 
demonstrates different power laws for densities of occupied (solid squares) and empty (light 
squares) states close to Fermi level and their common Gaussian asymptotics far from it. b) 
Single-particle density of states at T = and x = 1/2. 

of e(n, n') over all occupied n and empty n' in this configuration, and, if it is negative, perform 
the corresponding move. This procedure is repeated until such a configuration cm({ii}) is 
reached that is positive. Then cjn({n}) corresponds to a local equilibrium (with respect 
to single-particle moves) and the respective value of E c = E min [co({n})], a functional of the 
initial configuration, realizes a local minimum of energy. 

Next, the system is "shaken up", that is, for the same initial dopant configuration d({n}), 
a new arbitrary initial configuration c^({n}) is created, restricted again by the normalization 
condition, Eq. (||). Then it is found that, with the same equilibration process, some new 
local equilibrium CA/'({n}) and a new local minimum i^ in [co({n})] are obtained. The presence 
of various local minima is indicative of a non-ergodic structure of the phase space (with a 
discrete topology restricted to single-particle moves). Since there is a one-to-one correspondence 
between the initial state and the final state of local equilibrium, the whole phase space of the 
system gets divided into a number of attraction domains, each corresponding to a definite local 
minimum. The absolute energy minimum, corresponding to the "true" ground state, can be 
defined as: 

E^= min E mm [c ({n})]. (3) 

co({n}) 

The consequent "shake ups" and equilibrations define a sort of Monte-Carlo process to approach 
the ground state and this numeric process can be completed within a reasonable time for 
not too big system (this is evidenced by the uniqueness of the corresponding configuration 
c^ s ({n})). As a co-product, the process also provides the "spectrum" of local minima ui 0C (E) = 
(S(E — E min [co])) co , a new characteristics of the strongly interacting disordered system. At a 
given dopant configuration d({n}), the single-particle and the pair spectra, z/ s _ p (e) and f p (e), 
are calculated with respect to each local minimum, including the true ground state. Within 
accuracy to statistical noise, there is no difference found among the curves taken with respect to 
different minima. This shows that the "true" ground state is not any outstanding point between 
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Figure 3: Power law exponents for densities of occupied (filled circles) and empty (light circles) 
single-particle states as functions of concentration of dopants. 

other equilibrium points. Finally all the results are averaged over various dopant configurations 
d({n}). For numeric simulations of the system, Eqs. (|l|, ||), we began from cubic samples of 
increasing size L, then a size independent behavior corresponding to thermodynamic limit is 
reached already at L ~ 10, and these averaged characteristics are presented in Figs. 2-5. 

The single-particle spectra u s ^ p (e) at different dopings, shown in Fig. 2a, b, reveal a well- 
defined Coulomb gap around Fermi energy e F , while their asymptotics far from e F is well 
described by the Gaussian law: ^ s - p (e) oc exp [— a(e — e F ) 2 }, in agreement with the known 
results of Lifshitz's theory of optimal fluctuation for "tail" states in disordered systems [16]. 

A new notable feature of these spectra is a pronounced asymmetry between the densities of 
empty and occupied states near e F , which are described by different power laws: 

1 \ f ( e ~ £f) 2+v \ e > e F , A , 
u -*® x \ \ e -e F y-v\ e<, (4) 

Here the asymmetry factors 7/ and rj" are not universal, but vary with the doping and tend to 
zero at x — > 1/2 (Fig. |3|). The latter fact can be easily understood as a consequence of the 
symmetry between filled and empty sites at this concentration. 

Hence the mean-field quadratic law is found to be exact only at half-filling (though the 
non-universal corrections to it, due to the higher order correlations, can be really small in the 
case of traditional doped semiconductors with r a). 

Another qualitative difference from the mean-field behavior was found in that the density of 
pair excitations v p (e) does not remain constant but tends to zero with e — > (Fig. £|). Notably, 
this result is concentration independent. It is of interest for the future studies, to check also 
the low-energy behavior for the spectra of higher order excitations. 

At least, the "spectrum" of local energy minima for the doping level x = 1/3 is shown 
in Fig. |j. Close to the value of absolute energy minimum (which is proportional to 

the sample volume), this distribution exhibits quadratic energy dependence (dashed line) and, 
farther from E^? n , it falls down by the Gaussian law (dotted line): oc exp [—a(E — E^ n ) 2 ] , 
(while the overall distribution width is independent of the sample size and relatively small). 

Besides the simplest cubic shape of the samples, we also examined the slabs with L\ = L 2 = 
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Figure 4: (a) Density of pair excitations at two different concentrations of dopants, 
(b) Distribution of local minima of full electronic energy (at x = 1/3). The dashed line marks 
a parabolic dependence and the dash-dotted curve is Gaussian. 



6 and L 3 = 30. The essential features of single-particle spectrum for them (overall extension 
and gap asymmetry) at x = 1/3 were found practically identical to those for cubic samples. At 
least, a modification of the above slab configuration was considered, realizing a "topologically 
closed" Euclidean bar (Fig. |5|) where each vector (m, n 2) L + n 3 ) is identified with (ni,n.2,n 3 ). 
This form is appropriate for direct simulations of current flow in a closed circuit under external 
electric field applied along e 3 . However in this case a special care should be taken to assure the 
continuity of total energy at transitions of particles through the "topological" interface: L 1, 
which must be equivalent to any "normal" interface n 3 «-> n 3 + 1. To this end, the interaction 
potential u(\n — m|) in Eq. ([!]) should be replaced by the modified potential: 

u(n, m) = w(|n — m|) + u(\n + Le 3 — m|), (5) 

(n 3 < m 3 ), related to the two "distances" shown in Fig. [5|. Evidently, this modified interaction, 
Eq. (H), which connects each pair of sites in two possible ways, turns to the common Coulomb 
law in the limit L — > 00. Numeric simulations for a "closed" (6,6,30) slab with the modified 
interaction, Eq. (^|), showed no significant difference in the single-particle spectrum, compared 
to the similar "open" slab and usual interaction u(r). Thus, the characteristics of the ground 
state and excitation spectra found from our simulations are not sensitive to the sample size, 
shape and topology and should correspond to the true thermodynamic limit of the strongly 
interacting disordered system, Eq. ([!]). 



3 Finite temperatures 

At finite temperatures, the system dynamics is determined by the thermally activated hops 
between nearest neighbor sites. As usually, these hops are considered to be controlled by the 
electron-phonon interaction in approximation of deformation potential [|15|]. If we consider only 
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Figure 5: A lattice sample "topologically closed" in one dimension (labeled from 1 to L), here 
two "distances" are defined for each pair of points n and m. 

longitudinal phonons with Debye dispersion, the transition rate between the neighboring sites n 
and n' with the pair energy difference e(n, n' = e), including the probabilities for both processes 
with phonon absorption (e > 0) and emission (e < 0), is given by a simple expression: 

e - e D sin(e/e D ) 2 2 
7(6) = 70 exp(e/r) - 1 ° {eD ~ e) (6) 

where = Tis/ais the Debye energy (for sound velocity s) and 70 is some constant proportional 
to the small tunneling matrix element. The presence of the sine term in Eq. (jy) is due to the 
extremely short range of hops (by one lattice constant) and this factor essentially reduces the 
transition amplitude at low energies, compared to the usual case of electron-phonon interaction 
in metals. 

If the mean interval of time between two consecutive transitions, in some volume where the 
correlations are sensible, is much longer than the transition time r tr itself (this is reasonable for 
sufficiently small 70 and atomically fast T tr ), we can consider that only one transition occurs 
at a time. Also, accordingly to the analysis by Knotek and Pollak [[Uj, for so strictly localized 
states we can neglect transitions where more than one electron participate. However, the 
difficulty with applying Eq. (||) directly to our system consists in the fact that the transition 
energy for a given pair of sites (in a given transition channel) is not a priori defined but 
depends on the overall configuration and is varied at transitions between other sites. Since 
those transitions occur at random moments of time, the transition rate for any given channel is 
itself a random function of time, correlated with all other channels. In mathematical language, 
this means a realization of non-Markovian branching random process [|17[]. To manage this 
problem numerically, we used the following algorithm. 

For any initial configuration <i({n}), c({n}), Eq. (^) defines the transition rates 7 nn / = 
7[e(n, n')] for all appropriate pairs n — ► n', and these alternatives (the transition channels) can 
be considered independent. In this approach, only one of the channels can be chosen for each 
consecutive transition, and this choice is simulated by doing the independent statistical trials 
for random transition times r n n /, accordingly to the rates 7 n , n ', and choosing the shortest time 
from the trial outputs. 

In each particular trial, a random number £, < £ < 1, is generated, producing the 
associated random "phase" (f> = — ln£. It relates to the random output transition time r at a 
constant transition rate 7 through: r = 0/7. If £ is distributed uniformly: = 9(^)9(1 — £), the 
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Figure 6: Realization of a transition through a certain channel with a random apriori value of 
phase and transition rates 7o,...,7i which depend on time in a random way. The moment 
t[ corresponds to the virtual value for this channel obtained at Tj_i; in fact this value was not 
realized because it "lost" the competition to a shorter time 7j obtained at Tj_i for other channel. 
The true transition time for this channel, r, is defined by the "win" of the output (f>iji in the 
competition with outputs for all other channels at 7$. 

related distribution for the phase is: P^ = exp(— 0), and hence the distribution for transition 
times: P T = 7exp(— 77"). The relation between the random phase and transition time can 
be generalized to the case when the transition rate 7 is not constant but a posteriori definite 
function of current time 7 = j(t): 

<P= [ l(t)dt (7) 
Jo 

(in our specific case this function is stepwise, see Fig. ^J). Since the integrand in the r.h.s. of Eq. 
(0) is nonnegative, there always exists a single finite solution for the transition moment (Fig. 
^ . This value is determined not only by the trial phase value but also by all the intermediate 
times Ti and transition rates 7, = 7(r, < r < r i+ i) expressing complicate correlations between 
the given transition and the preceding ones. Each time interval Arj = r i+ i — Ti is determined 
by the result of competition between the virtual values for all the channels j possible after the 
moment 77 

An = min(0 J °' ) /7?' ) ) (8) 
j 

(i) 

where 7 4 is the transition rate in j-th channel between the moments Tj and Tj+i, and: 

(4) = <j>® - [ 1 >yV\t)dt (9) 

is the "residual phase at T{ for j-th channel. This channel is opened at a certain time moment 
t^\ when it is given the initial random phase value according to P^. Then is con- 
secutively reduced to 4\ \ by Eq. (§), at each intermediate transition, until such a moment 
Ti is reached that the virtual value Mi f° r this channel "wins" the competition, Eq. (§). 
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Figure 7: Temporal evolution of full energy in a cubic sample with L = 8 and x = 1/3 at 
different temperatures. Time units are cdIq 1 x 10 5 and energy units are e 2 /na. Inset: relaxation 
time t r and statistical entropy S, Eq. 11, vs temperature; a rapid increase of t r below T ~ 0.05 
indicates the freezing of a glassy system. 

(i) (i) 

Then the residual phase value attained at Tj+i = r, + 4>i n{ is just zero, corresponding to 
an exact solution of Eq. (0) at r = r i+1 . After a particle has performed a transition through 
j-th channel, this channel (together with the whole set of channels j, having common initial 
site with j) is considered closed, and a number of new channels is opened for all empty sites, 
neighbors to the new occupied site. Also, the change of the system configuration produced by 
the transition in j-th channel implies all the rates 7,- , j' ^ j, to be changed for some new 
values 7^1, and, after opening of new channels with corresponding initial phases and transition 
rates, the whole process is continued. At each z-th transition in the system, the values of the 
energy transfer (either positive or negative) and of the time interval A-Tj, past the preceding 
transition are recorded. 

It is important that, unlike a system of independent carriers, the whole passage from opening 
to closure of a channel is not typically a pair process and it is related to the overall density of 
many-body states. 

Examples of time records for the full energy as functions of full time r» = Y^i>=\ ^ r «' j a t 
different values of temperature T and concentration x = 1/3 for a cubic sample with L — 8, 
are shown in Fig. [7|. They all demonstrate an initial rapid descent to the regime of dynamical 
equilibrium, within a certain relaxation time t r . This time is almost independent of temperature 
except at very low temperatures (below ~ 0.05, in our energy units e 2 /na) when t r grows very 
rapidly (see inset to Fig. |7|). The latter can serve as an indication of the "freezing" process 
in the glassy system, though, of course, there is no sharply defined critical temperature in a 
finite size sample. In the equilibrium regime, both the mean energy (E) and its dispersion 
5E = ((E 2 ) — (E) 2 ) 1 / 2 are obtained as certain functions of temperature. 

Since the volume of our system is kept fixed and the average in time of total energy (E) as a 
function of temperature is known, the specific heat C v can be readily obtained by differentiation: 
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Figure 8: Specific heat vs temperature for the system corresponding to Fig. 8. Points stand 
for the numerically calculated derivative dE/dT, and the solid line is a guide for the eye. 



C v = d(E)/dT ||. 

The corresponding numeric result is shown in Fig. [8|. It displays a linear behavior at low 
temperatures: C v oc T, characteristic both of the Fermi-liquid systems p(J and of the glassy 



systems [21 



Another distinctive feature of the considered specific heat is a pronounced maximum at 
T ~ 2e/), corresponding to saturation of the relaxation channels above the Debye temperature. 
Through the relation: 

S(T) = [ T~ l C v {T)dT (10) 

Jo 

the thermodynamical entropy S(T) is fully determined by the function (E(T)}. On the other 
hand, S is also related to 5E pOj] : 

S = ln[u m _ b ((E))5E], (11) 

which permits to estimate the total density of many-body states v m _b(E), very difficult in other 
approaches. From the comparison of Eqs. ([TOD and ([TTD at T — > 0, we conclude that u m -b(E) 
attains a finite value near E^ n , though, for a quantitative accuracy, one perhaps needs more 
precise and detailed data on 5E(T) and (E(T)) than those in Figs. [7|, [|. 

Next, the temperature behavior of the hopping conductivity cr(T), accordingly to the Ein- 
stein relation: a = ne 2 D /(fc B T) ^2|, can be estimated from that of the diffusion coefficient 
D. Since all the hops have a standard length a, the diffusion coefficient D = a 2 /(3r ) is fully 
determined by the mean lifetime To of a localized state, and the latter is merely the inverse 
of the average number of hops per one particle per unit time: r = iV limj_, 00 (rj/i). Then, 
from the plots of full number of hops i vs full time T{ at different temperatures and the same 
concentration x = 1/3 (they all turn perfectly linear after the same relaxation time t r as that 
for energy, Fig. 9a), we deduced the values of D, proportional to the slopes di/dri. Finally, the 
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Figure 9: a) The slopes of linear functions "full number of hops vs full time" provide the 
temperature dependence of the diffusion coefficient, b) Double logarithmic plot for conductivity 
a vs In (1/T) at x — 1/3. The slope is different from the mean-field value 1/2. 

double logarithmic plot: ln(ln cr — In a) vs ln(l/T), Fig. 9b) (where the fitting parameter cr 
was adjusted to get the best linearity), permits to infer the modified hopping conductivity law: 

a(T) = a o exp[-(T /T)% (12) 

with a ~ 0.87 and reasonably low To ~ 0.1. It is of interest to compare these figures with 
the mean-field values ||: a = 1/2 and T ~ 2.1 (the latter is obtained using the estimate for 
localization radius r ~ 0.4a by the Mott criterion at x = 1/3, see Introduction). 

4 Conclusions 

A model of strongly disordered lattice system with long-range Coulomb interactions between 
localized charge carriers has been considered. The total electronic energy is characterized by 
the presence of multiple metastable minima (including the true ground state), and different 
types of excitation spectra over these minima. A numeric procedure, accounting for all many- 
body correlations in finite size samples, confirms the existence of Coulomb gap in the single- 
particle spectrum and also provides corrections to the known mean-field theory results, as 
asymmetry of the gap at non-half-filling, vanishing density of low energy pair excitations, 
modified temperature exponent for hopping conductivity. The further analysis of this model 
can involve direct simulations of current flow in a "topologically closed" sample (Sec. 2) and 
formation of cluster states at finite tunneling amplitudes between nearest neighbour sites with 
sufficiently small pair energies. 

This research was supported by the Portuguese program PRAXIS XXI through the project 
2/2.1/FIS/302/94 and under personal Grants BPD 14226/97 (Yu.G.P.) and BM/12717/97 (J. 
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